### Setup
library(hash)
rm(list=ls())
source("reg-pic.R")
source("postpred-pics.R")
load("../fits/main.RData")

## A little function to calculate the mode
dmode <- function (x) density(x)$x[which.max(density(x)$y)]

## Setup posterior pred case
## Use means for every variable except those listed as keys.  We used
## median ingov and look at big 3 groups, and use German candidates as
## our exemplars

case.hash=hash(
    "ingov" = median,
    "big3" = function (x) 1,
    "Germany" = function (x) 1,
    "Czech Republic" = function (x) 0,
    "France" = function (x) 0,
    "Greece" = function (x) 0,
    "Hungary" = function (x) 0,
##    "Ireland" = function (x) 0,
    "Italy" = function (x) 0,
    "Netherlands" = function (x) 0,
    "Romania" = function (x) 0,
    "Spain" = function (x) 0,
    "UK" = function (x) 0)

### Regression plot

plotplist(post[[1]],
          vars=c("(int)", "emph.europe", "big3", "ingov", "propseats",
                  "pro_anti_eu"),
          vars.names=c("Intercept", "EU Emphasis", "Big 3", "In Govt.",
              "% Seats", "Pro/Anti EU"),
          base.group="inc",
          other.groups=c("non", "reg", "nat"),
          group.names=c("0", "R", "N"))


### Posterior prediction plotsxs

## EU Emphasis
plot.plist.postpred.contin(post[[1]], diag(4), "emph.europe",
                           c("No Experience",
                             "Local/Regional",
                             "National",
                             "Incumbent"),
                           case.hash,
                           xlab="Emphasis on Europe",
                           point.col=c("red", "cyan", "blue", "orange"),
                           hpd.col=c(rgb(255, 0, 0, 50, maxColorValue=255),
                                     rgb(0, 255, 255, 50, maxColorValue=255),
                                     rgb(0, 0, 255, 50, maxColorValue=255),
                                     rgb(255, 133, 0, 50, maxColorValue=255)),
                           hpd.thin=3, hpd.prob=0.95, point.func=mean,
                           legend.x=60, legend.y=0.55, legend.horiz=F,
                           legend.bty='n', point.lwd=4,
                           axis.labels=round(seq(min(y$emph.europe),
                                           max(y$emph.europe),
                                           length.out=6)))

## Big3
p<-plot.plist.postpred.dummy(post[[1]], diag(4), "big3",
                          c("No Experience",
                            "Local/Regional",
                            "National",
                            "Incumbent"),
                          case.hash,
                          xlab="Group Membership",
                          point.col=c("red", "cyan", "blue", "orange"),
                          point.lwd=4,
                          hpd.col=c("red", "cyan", "blue", "orange"),
                          hpd.prob=0.95,
                          point.func=mean,
                          legend.bty='n',
                          axis.labels=c("Other Group", "Big 3"),
                          legend.x=-.5, legend.y=1, legend.horiz=F,
                          return.probs=T)

## Seats in national leg
plot.plist.postpred.contin(post[[1]], diag(4), "propseats",
                           c("No Experience",
                             "Local/Regional",
                             "National",
                             "Incumbent"),
                           case.hash,
                           xlab="% Seats in National Legislature",
                           point.col=c("red", "cyan", "blue", "orange"),
                           hpd.col=c(rgb(255, 0, 0, 50, maxColorValue=255),
                                     rgb(0, 255, 255, 50, maxColorValue=255),
                                     rgb(0, 0, 255, 50, maxColorValue=255),
                                     rgb(255, 133, 0, 50, maxColorValue=255)),
                           hpd.thin=3, hpd.prob=0.95, point.func=mean,
                           legend.bty='n', point.lwd=4,
                           legend.x=70, legend.y=1, legend.horiz=F)


### Various statistics quoted in the text

## Calc prob coefficients less than 0 for big3
sum(contrast(post[[1]], "big3", "nat", "inc")<0)/5000

## Prob change for selecting inc
summary(p[[4]][,2]-p[[4]][,1])
HPDinterval(as.mcmc(p[[4]][,2]-p[[4]][,1]))

## Prob big3 selects inc over nat
summary(p[[4]][,2]-p[[3]][,2])

## Prob non-big3 selects inc over nat
summary(p[[4]][,1]-p[[3]][,1])

## Prob big3 premium on incumbency bigger than non-big3
sum((p[[4]][,2]-p[[3]][,2]) > (p[[4]][,1]-p[[3]][,1]))/5000

## Prob ingov coefs less than zero
sum(contrast(post[[1]], "ingov", "non", "inc")<0)/5000
sum(contrast(post[[1]], "ingov", "reg", "inc")<0)/5000
sum(contrast(post[[1]], "ingov", "nat", "inc")<0)/5000